arquivos <- list.files("Dados", pattern = "*.csv", full.names = TRUE)

dados <- list()
for (i in seq_along(arquivos)) {
  dados[[i]] <- read.csv(arquivos[i], sep=';')
}

dados_combinados <- do.call(rbind, dados)

dados_selecionados = subset(dados_combinados,
                            dados_combinados$Municipio == 'GOIANIA' &
                              dados_combinados$Produto == 'GASOLINA' &
                              dados_combinados$Revenda == 'POSTO XODO LTDA')

dados_selecionados <- dados_selecionados[, c(4, 11, 12, 13)]
head(dados_selecionados)
##               Revenda  Produto Data.da.Coleta Valor.de.Venda
## 1365  POSTO XODO LTDA GASOLINA     03/01/2018          4,499
## 1368  POSTO XODO LTDA GASOLINA     03/01/2018          4,499
## 19392 POSTO XODO LTDA GASOLINA     10/01/2018          4,599
## 19395 POSTO XODO LTDA GASOLINA     10/01/2018          4,499
## 38348 POSTO XODO LTDA GASOLINA     17/01/2018          4,599
## 38351 POSTO XODO LTDA GASOLINA     17/01/2018          4,499
#table(dados_selecionados$Revenda)

dados3 <- dados_selecionados

rm(dados_combinados, dados)

1 Resumo

  • Período: 2018 a 2023;
  • Posto X;
  • Quantidade de observações iniciais: 365
  • Método de agrupamento: média mensal;
  • Quantidade de observações finais: 72;
  • Imputação de dados na observação dos meses listados abaixo foi feita a partir da média daquele ano;
  • Tentativa de ajuste com um modelo SARIMA;

A partir da base de dados dos preços de combustíveis na cidade de Goiânia, será feita a filtragem listada acima. Foi selecionado um único posto revendedor - o posto com maior número de observações no período citado - contendo 365 valores divididos em 188 datas distintas; a seguir, foi tomada a média de cada mês essa será a variável estudada; foi feito a imputação de valores a cinco meses que não possuíam dados (nov-2018, dez-2018, fev-2019, set-2020, fev-2023) por meio da média dos seus respectivos anos; por fim,

2 Organização dos dados

2.1 Antes

dados = dados3g[, -c(1,2)]

#head(dados)

summary(dados)
##  Data.da.Coleta       Valor.de.Venda      mes                ano           
##  Min.   :2018-01-03   Min.   :3.699   Length:365         Length:365        
##  1st Qu.:2019-12-04   1st Qu.:4.690   Class :character   Class :character  
##  Median :2021-04-07   Median :4.990   Mode  :character   Mode  :character  
##  Mean   :2021-03-13   Mean   :5.396                                        
##  3rd Qu.:2022-09-29   3rd Qu.:5.890                                        
##  Max.   :2023-12-26   Max.   :7.870

Nota-se que o comportamento dos preços do posto X se assemelha ao geral.

2.2 Depois

# Meses faltantes
m <- check_missing_months(dados, "Data.da.Coleta")
## [1] 5
print(m)
## [1] "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# xodo: 5 faltas "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# hiper moreira: 18 faltas
# viena: 5
# Monaco: 4

Como há 5 meses sem nenhuma observação foi utilizado da imputação de dados. Para cada mês listado acima será imputado o valor da média do respectivo ano, portanto, novembro e dezembro de 2018 terão o mesmo valor.

# Calculando médias anuais
aux = dados %>% group_by(ano) %>% mutate(media_anual = mean(Valor.de.Venda))

#unique(aux[, c(4, 5)])

# "2018-11" "2018-12" "2019-02" "2020-09" "2023-02"
# Adicionar os seguintes valores
linha1 = data.frame(matrix(c("2018-11-01", 4.661141, "2018-11", "2018",
                             "2018-12-01", 4.661141, "2018-12", "2018",
                             "2019-02-01", 4.665063, "2019-02", "2019",
                             "2020-09-01", 4.460714, "2020-09", "2020",
                             "2023-02-01", 5.702877, "2023-02", "2023"),
                           byrow= T, ncol=4))
colnames(linha1) <- colnames(dados)

linha1[, 1] <- as.Date.character(linha1[, 1])
linha1[, 2] <- as.numeric(linha1[, 2])
d1 = unique(dados_imput[, c(3, 5)])

qcc::qcc(d1$media_mensal, type = 'xbar.one',
         title = 'Carta da média amostral (Dados finais)')

## List of 11
##  $ call      : language qcc::qcc(data = d1$media_mensal, type = "xbar.one", title = "Carta da média amostral (Dados finais)")
##  $ type      : chr "xbar.one"
##  $ data.name : chr "d1$media_mensal"
##  $ data      : num [1:72, 1] 4.54 4.52 4.5 4.4 4.38 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:72] 4.54 4.52 4.5 4.4 4.38 ...
##   ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
##  $ sizes     : int [1:72] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center    : num 5.34
##  $ std.dev   : num 0.19
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] 4.76 5.91
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"
qcc::ewma(d1$media_mensal, title = 'Carta MMEP (Dados finais)')

## List of 15
##  $ call      : language qcc::ewma(data = d1$media_mensal, title = "Carta MMEP (Dados finais)")
##  $ type      : chr "ewma"
##  $ data.name : chr "d1$media_mensal"
##  $ data      : num [1:72, 1] 4.54 4.52 4.5 4.4 4.38 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:72] 4.54 4.52 4.5 4.4 4.38 ...
##   ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
##  $ sizes     : int [1:72] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center    : num 5.34
##  $ std.dev   : num 0.19
##  $ x         : int [1:72] 1 2 3 4 5 6 7 8 9 10 ...
##  $ y         : Named num [1:72] 5.18 5.05 4.94 4.83 4.74 ...
##   ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
##  $ sigma     : num [1:72] 0.0381 0.0487 0.0545 0.0579 0.0599 ...
##  $ lambda    : num 0.2
##  $ nsigmas   : num 3
##  $ limits    : num [1:72, 1:2] 5.22 5.19 5.17 5.16 5.16 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ violations: Named int [1:70] 1 2 3 4 5 6 7 8 9 10 ...
##   ..- attr(*, "names")= chr [1:70] "1" "2" "3" "4" ...
##  - attr(*, "class")= chr "ewma.qcc"
qcc::cusum(d1$media_mensal, title = 'Carta "Soma cumulativa" (Dados finais)')

## List of 14
##  $ call             : language qcc::cusum(data = d1$media_mensal, title = "Carta \"Soma cumulativa\" (Dados finais)")
##  $ type             : chr "cusum"
##  $ data.name        : chr "d1$media_mensal"
##  $ data             : num [1:72, 1] 4.54 4.52 4.5 4.4 4.38 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics       : Named num [1:72] 4.54 4.52 4.5 4.4 4.38 ...
##   ..- attr(*, "names")= chr [1:72] "1" "2" "3" "4" ...
##  $ sizes            : int [1:72] 1 1 1 1 1 1 1 1 1 1 ...
##  $ center           : num 5.34
##  $ std.dev          : num 0.19
##  $ pos              : num [1:72] 0 0 0 0 0 0 0 0 0 0 ...
##  $ neg              : num [1:72] -3.66 -7.43 -11.33 -15.75 -20.24 ...
##  $ head.start       : num 0
##  $ decision.interval: num 5
##  $ se.shift         : num 1
##  $ violations       :List of 2
##  - attr(*, "class")= chr "cusum.qcc"
# AGRUPAMENTO TRIMESTRAL
m1 = matrix(d1$media_mensal, ncol = 3, byrow = T)
m11 = matrix(d1$mes, ncol=3, byrow = T)

qcc::qcc(data = m1, type = "S")

## List of 11
##  $ call      : language qcc::qcc(data = m1, type = "S")
##  $ type      : chr "S"
##  $ data.name : chr "m1"
##  $ data      : num [1:24, 1:3] 4.54 4.4 4.8 4.98 4.57 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:24] 0.0223 0.2833 0.1042 0.1835 0.0876 ...
##   ..- attr(*, "names")= chr [1:24] "1" "2" "3" "4" ...
##  $ sizes     : int [1:24] 3 3 3 3 3 3 3 3 3 3 ...
##  $ center    : num 0.184
##  $ std.dev   : num 0.207
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] 0 0.471
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"
qcc::qcc(data = m1, type = "R")

## List of 11
##  $ call      : language qcc::qcc(data = m1, type = "R")
##  $ type      : chr "R"
##  $ data.name : chr "m1"
##  $ data      : num [1:24, 1:3] 4.54 4.4 4.8 4.98 4.57 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:24] 0.0444 0.4976 0.1809 0.3179 0.1751 ...
##   ..- attr(*, "names")= chr [1:24] "1" "2" "3" "4" ...
##  $ sizes     : int [1:24] 3 3 3 3 3 3 3 3 3 3 ...
##  $ center    : num 0.341
##  $ std.dev   : num 0.202
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] 0 0.878
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"

A aplicação das cartas aos dados autocorrelacionados indica uma quantidade muito grande de pontos fora de controle. É de se esperar que o processo não seja bem medido e tais dados façam com que as cartas tornem-se, frequentemente, errôneas.

3 Série temporal

par(bg='gray')
d1 = unique(dados_imput[, c(3, 5)])

serie = ts(d1[, 2], start=c(2018, 1), frequency = 12)
#serie = ts(dados_imput$media_mensal, start=c(2018, 1), end=c(2023, 12), frequency = 12)
plot(serie, ylab="Valor", xlab="Obs.")

forecast::autoplot(decompose(serie))

Verificando a presença de autocorrelação entre as observações da série nota-se que há o decaimento exponencial à medida que se aumenta o lag.

par(bg='gray')
acf(d1$media_mensal, main = "Correlação entre as médias mensais")

4 Modelo ARIMA

Após toda organização do conjunto de dados, foram testados diversos modelos SARIMA. A seguir há uma relação dos modelos que apresentaram resíduos com: normalidade, estacionariedade e independência.

Outros modelos que também apresentaram todos pressupostos.

#aux = forecast::auto.arima(serie)

ar_ = forecast::Arima(serie, order = c(2, 3, 0),
                      seasonal = c(0,0,0), method = "ML")

summary(ar_)
## Series: serie 
## ARIMA(2,3,0) 
## 
## Coefficients:
##           ar1      ar2
##       -0.9819  -0.4259
## s.e.   0.1076   0.1073
## 
## sigma^2 = 0.2455:  log likelihood = -48.96
## AIC=103.93   AICc=104.3   BIC=110.63
## 
## Training set error measures:
##                       ME      RMSE      MAE       MPE    MAPE      MASE
## Training set 0.001505884 0.4779935 0.353291 0.1276067 6.67558 0.3425822
##                     ACF1
## Training set -0.09204694
# AMBOS PARÂMETROS DO MODELO ARIMA SÃO NEGATIVOS!!

shapiro.test(residuals(ar_))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(ar_)
## W = 0.96813, p-value = 0.06461
tseries::adf.test(residuals(ar_))
## Warning in tseries::adf.test(residuals(ar_)): p-value smaller than printed
## p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  residuals(ar_)
## Dickey-Fuller = -6.1038, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
Box.test(residuals(ar_), type="Ljung")
## 
##  Box-Ljung test
## 
## data:  residuals(ar_)
## X-squared = 0.63581, df = 1, p-value = 0.4252
forecast::checkresiduals(ar_)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,3,0)
## Q* = 23.918, df = 12, p-value = 0.02087
## 
## Model df: 2.   Total lags used: 14

5 Aplicação das cartas de controle

#agrup = qcc::qcc.groups(as.vector(ar_$residuals),
#                        as.vector(unique(dados_imput[, 3])))

g1 = qcc::qcc(ar_$residuals, type = "xbar.one", plot = T,
              xlab="Grupo (Mês)",
              ylab="Estatística do grupo",
              add.stats = F,
              label.limits = c("LCI", "LCS"))

g1$violations$beyond.limits
## [1] 57
g2 = qcc::ewma(ar_$residuals)

g3 = qcc::cusum(ar_$residuals)

#plot(g1, restore.par = F)
#head(d1)
d1[order(d1$media_mensal, decreasing = T), ]
## # A tibble: 72 × 2
## # Groups:   mes [72]
##    mes     media_mensal
##    <chr>          <dbl>
##  1 2022-05         7.63
##  2 2022-04         7.62
##  3 2022-06         7.61
##  4 2022-03         7.53
##  5 2022-01         7.36
##  6 2021-11         7.34
##  7 2021-12         7.03
##  8 2021-10         7.02
##  9 2022-02         6.99
## 10 2021-09         6.66
## # ℹ 62 more rows
d1[56:58, ]
## # A tibble: 3 × 2
## # Groups:   mes [3]
##   mes     media_mensal
##   <chr>          <dbl>
## 1 2022-08         5.10
## 2 2022-09         4.96
## 3 2022-10         4.99
g1$statistics[56:58]
##        56        57        58 
## 0.9355005 1.6073610 0.3819759
g1 = qcc::qcc(ar_$residuals, type = "xbar.one", plot = F,
              add.stats = F)

# 1. Convert qcc object to a data frame for ggplot2
qcc_data <- data.frame(
  x = 1:length(g1$data), # x-axis (sample number)
  Value = g1$data,      # y-axis (data values)
  CL = mean(g1$data[, 1]),     # Center Line
  UCL = g1$limits[2],    # Upper Control Limit
  LCL = g1$limits[1],    # Lower Control Limit
  violations = ifelse(g1$data > g1$limits[2] | g1$data < (g1$limits[1] - 3*(g1$limits[2] - g1$limits[1])), "Violação", "Sob controle"))

# 2. Create the ggplot2 plot
gg1 <- ggplot(qcc_data, aes(x = x, y = Value)) +
  theme_classic() +
  geom_point(aes(color = violations)) +  # Points, colored by violations
  geom_line() +                      # Connect the points
  geom_hline(yintercept = qcc_data$CL[1], linetype = "solid", color = "black",
             linewidth = .1) + # Center line
  geom_hline(yintercept = qcc_data$UCL[1], linetype = "dashed", color = "red",
             linewidth = .3) + # Upper control limit
  geom_hline(yintercept = qcc_data$LCL[1], linetype = "dashed", color = "blue",
             linewidth = .3) + # Lower control limit
  labs(title = "", x = "Grupo (Mês)", y = "Estatística do grupo") + # Labels
  scale_color_manual(values = c("Violação" = "red", "Sob controle" = "black")) + # Custom colors
  annotate(geom = 'text', label = c(paste0("LSC = ", round(qcc_data$UCL[1], 3)),
                                    paste0("LC = ", round(qcc_data$CL[1], 3)),
                                    paste0("LIC = ", round(qcc_data$LCL[1], 3))),
           size = 3, x = c(5,1,5), y = c(qcc_data$UCL[1] + .1 ,
                                         qcc_data$CL[1] + .1,
                                         qcc_data$LCL[1] + .1)) +
  theme(legend.position = "none") # Remove the legend

plotly::ggplotly(gg1)

5.1 Adicionais

ar1 = matrix(ar_$residuals, ncol = 3, byrow = T)

ar11 = matrix(d1$mes, ncol = 3, byrow = T)
ar111 = matrix(c(1:72), ncol = 3, byrow = T)

am1 = qcc::qcc(ar1, type = "R")

g1$violations$beyond.limits
## [1] 57
#ar11[g1$violations$beyond.limits, ]

#qcc::qcc(ar1, type="S")
ar1 = matrix(ar_$residuals, ncol = 4, byrow = T)

ar11 = matrix(d1$mes, ncol = 4, byrow = T)
ar111 = matrix(c(1:72), ncol = 4, byrow = T)

am1 = qcc::qcc(ar1, type = "R")

am1$violations$beyond.limits
## [1] 14 15
#ar11[am1$violations$beyond.limits, ]

qcc::qcc(ar1, type="S")

## List of 11
##  $ call      : language qcc::qcc(data = ar1, type = "S")
##  $ type      : chr "S"
##  $ data.name : chr "ar1"
##  $ data      : num [1:18, 1:4] 0.00104 0.10152 0.28529 0.16085 0.04526 ...
##   ..- attr(*, "dimnames")=List of 2
##  $ statistics: Named num [1:18] 0.0244 0.4928 0.363 0.2921 0.2695 ...
##   ..- attr(*, "names")= chr [1:18] "1" "2" "3" "4" ...
##  $ sizes     : int [1:18] 4 4 4 4 4 4 4 4 4 4 ...
##  $ center    : num 0.457
##  $ std.dev   : num 0.496
##  $ nsigmas   : num 3
##  $ limits    : num [1, 1:2] 0 1.04
##   ..- attr(*, "dimnames")=List of 2
##  $ violations:List of 2
##  - attr(*, "class")= chr "qcc"
ar1 = matrix(ar_$residuals, ncol = 3, byrow = T)

ar11 = matrix(d1$mes, ncol = 3, byrow = T)
ar111 = matrix(c(1:72), ncol = 3, byrow = T)

am1 = qcc::qcc(ar1, type = "S")

g1$violations$beyond.limits
## [1] 57
#ar11[g1$violations$beyond.limits, ]
# 1. Convert qcc object to a data frame for ggplot2
qcc_data <- data.frame(
  x = 1:length(am1$statistics), # x-axis (sample number)
  Value = am1$statistics,      # y-axis (data values)
  CL = mean(am1$statistics),     # Center Line
  UCL = am1$limits[2],    # Upper Control Limit
  LCL = am1$limits[1],    # Lower Control Limit
  violations = ifelse(am1$statistics > am1$limits[2] | am1$statistics < (am1$limits[1] - 3*(am1$limits[2] - am1$limits[1])), "Violação", "Sob controle"))

# 2. Create the ggplot2 plot
gg1 <- ggplot(qcc_data, aes(x = x, y = Value)) +
  theme_classic() +
  geom_point(aes(color = violations)) +  # Points, colored by violations
  geom_line() +                      # Connect the points
  geom_hline(yintercept = qcc_data$CL[1], linetype = "solid", color = "black",
             linewidth = .1) + # Center line
  geom_hline(yintercept = qcc_data$UCL[1], linetype = "dashed", color = "red",
             linewidth = .3) + # Upper control limit
  geom_hline(yintercept = qcc_data$LCL[1], linetype = "dashed", color = "blue",
             linewidth = .3) + # Lower control limit
  labs(title = "", x = "Grupo (Mês)", y = "Estatística do grupo") + # Labels
  scale_color_manual(values = c("Violação" = "red", "Sob controle" = "black")) + # Custom colors
  annotate(geom = 'text', label = c(paste0("LSC = ", round(qcc_data$UCL[1], 3)),
                                    paste0("LC = ", round(qcc_data$CL[1], 3)),
                                    paste0("LIC = ", round(qcc_data$LCL[1], 3))),
           size = 3, x = c(5,1,5), y = c(qcc_data$UCL[1] + .1 ,
                                         qcc_data$CL[1] + .1,
                                         qcc_data$LCL[1] + .1)) +
  theme(legend.position = "none") # Remove the legend

plotly::ggplotly(gg1)